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THE SIMPLEST MIXED FINITE ELEMENT METHOD FOR LINEAR ELASTICITY 
IN THE SYMMETRIC FORMULATION ON n-RECTANGULAR GRIDS 
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» I , Abstract. A family of mixed finite elements is proposed for solving the first order system of linear elasticity 

^1^' equations in any space dimension, where the stress field is approximated by symmetric finite element tensors. 

■^^ ' This family of elements has a perfect matching between the stress components and the displacement. The 

discrete spaces for the normal stress era, the shear stress aij and the displacement Ui are span{l,Xi}, 
span{l, Xi, Xj} and span{l}, respectively, on rectangular grids. In particular, the definition remains the 
same for all space dimensions. As a result of these choices, the theoretical analysis is independent of the 
spatial dimension as well. In ID, this element is nothing else but the ID Raviart- Thomas element, which 
.^^ ' is the only conforming element in this family. In 2D and higher dimensions, they are new elements but of 

}—^ , the minimal degrees of freedom. The total degrees of freedom per element is 2 plus 1 in ID, 7 plus 2 in 2D, 

^-H i and 15 plus 3 in 3D. The previous record of the least degrees of freedom is, 13 plus 4 in 2D, and 54 plus 12 

r^ ' in 3D, on the rectangular grid. These elements are the simplest element for any space dimension. 

The well-posedness condition and the optimal a priori error estimate of the family of finite elements are 
proved for both pure displacement and traction problems. Numerical tests in 2D and 3D are presented to 
show a superiority of the new element over others, as a superconvergence is surprisingly exhibited. 
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o 

en ' 1- Introduction 

The first order system of equations, for tlie symmetric stress field ct G E := i/(div,i7,§) and the dis- 
placement field ueV :— L^(fi,]R"), reads: Find (ct, m) e S x y such that 

{Act, t) + (divr, u) = Vr e S, 

^.' ^^'^' (divfj, i;) = (/, w) \tv€V. 

Here the symmetric tensor-valued stress space S and the vector-valued displacement space V are, respec- 
tively, 

(1.2) i/(div, 17, §) = { {cj,,)^^^ e H{d\y, n) I a,, = a,. }, 

(1.3) L2(f],R") = {(«!, ■■■ ,u^f \u,eL^{n)y 

In ID, one example of the problem (ll.ip is the mixed formulation of the ID Poisson equation; In 2D and 
3D, the stress-displacement formulation based on the Hellinger-Reissner principle for the linear elasticity 
can be regarded as a celebrated example of (|l.ip . 

Because of the symmetry constraint on the stress tensor, (Jij = (jji, it is extremely difficult to construct 
stable conforming finite elements of (11.11) even if for 2D and 3D, as stated in the plenary presentation 
to the 2002 International Congress of Mathematicians by D. N. Arnold. Hence compromised works use 
composite elements [6l [22j, or enforce the symmetry condition weakly [2 [5l [TTJ [24j [27l [28l [29]. The 
landmarks in this direction are the respective works of Arnold and Winther [8 and Arnold, Awanou, and 



The first author was supported by the NSFC Project 11271035, and in part by the NSFC Key Project 11031006. 

1 



Vk 



JUN HU, HONGYING MAN, AND SHANGYOU ZHANG 





{l,x,y,y^,xy} 








{I,a:,y,x2,y2} 






K ■ 


{i,y} 



{l,x,y,x ,y } 



{l,x,y,xy,x'^} 



{1,4 



{1,4 



{l,a;,y} 



{1} 



{^,x,y} 



{hy} 



{1} 



Figure 1. 2D elements by Hu-Shi [^T], Yi [21] and this paper. 



Winther ^. In particular, a sufficient condition of the discrete stable method is proposed in these two 
papers, which states that a discrete exact sequence guarantees the stability of the mixed method. Based 
on such a condition, conforming mixed finite elements on the simplicial and rectangular triangulations are 
developed for both 2D and 3D [H |3l IH [8] . In order to keep conformity the vertex degrees of freedom are in 
particular employed in these conforming methods. To avoid the complexity of conforming mixed element 
and also vertex degrees of freedom, new weak-symmetry finite elements [71 [161 UHl US], non-conforming 
finite elements [S] [211 [13 [131 [201 IS] are constructed. See also (THl [TO] for the enrichment of nonconforming 
elements of [211 [23] to conforming elements. However, most of these elements are difficult to be implemented; 
numerical implementation can only be found in [131 [Ml [31] so far, all in 2D. 
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Figure 2. The 3D element of this paper. 



In this paper, a new family of minimal, any space-dimensional, symmetric, nonconforming mixed finite 
elements for the problem (|l.ip is constructed. It is motivated by a simple fact that, by p.2p . the derivative on 
a normal stress component (Th is only in Xi direction; while those on ct^j are only in Xi and Xj directions. Thus, 
the minimal finite element space for un would be span{l,a;i} on each n-dimensional rectangular element; 
the minimal finite element space for aij would be span{l, Xi, Xj} on each n-dimensional rectangular element. 
For the displacement (jl.31) , there is no derivative and the minimal finite element space would be the constant 
space span{l}. The spaces are displayed in the right diagram in Figure [Hand in Figure [2l Surprisingly, it is 
shown that these minimal finite element spaces can actually form a family of stable and convergent methods 
for (|l.ip . However, the analysis herein has to overcome the difficulty to prove the discrete inf-sup condition, 
one key ingredient for the stability analysis of the mixed finite element method 12 , and the difficulty 
related to nonconformity of the discrete spaces for the stresses. For both the elasticity problem and the 
Poisson problem, the stability analysis of mixed finite element methods in literature is established by special 
commuting properties of canonical interpolation operators defined by degrees of freedom of discrete stress 
spaces, see, for instance, [T1[31[31[S] and [H]. To overcome the first difficulty, a new macro-element technique 
is proposed to prove a Fortin Lemma for mixed methods under consideration. Note that the macro-element 
technique is widely used to analyze the stability of mixed methods for the Stokes problem, see [12] and 
references therein. However, it is not used to the elasticity problem before. For the pure displacement 
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problem, an explicit constructive proof is also given for the discrete inf-sup condition. In order to deal with 
the second difficulty, a superconvergence property of the consistency error is proved. The mathematical 
elegance and beauty of this family of minimal elements is gestated within, besides the perfect matching, 
the independence of the spatial dimension n. In n dimension, the constructive proof of the discrete inf-sup 
condition can be divided into n steps of that for the ID Raviart-Thomas element, and the consistency error 
can be decomposed as n two-dimensional consistency errors (For ID, there is no consistency error.) 

The superiority of the family of elements over the existing elements in the literature is its simplicity and 
high accuracy. In fact, a family of 2D rectangular, conforming elements, of which the lowest order has 45 
stress and 12 displacement degrees of freedom per element, is proposed in [^. A nonconforming mixed 
finite element based on rectangular grids is proposed with 19 stress and 6 displacement degrees of freedom 
on each element in |30j . Later on, a simplified mixed finite element on 2D rectangular grids is constructed 
with 13 stress and 4 displacement degrees of freedom on each rectangle independently in [2TJ |3T], see the 
left diagram in Figure [TJ which is the simplest rectangular element of first order in 2D in the literature so 
far. Doubtless, the 2D element with 7 stress and 2 displacement degrees of freedom on each rectangle of 
this paper is the simplest rectangular element, see the right diagram in Figure [T] Due to a perfect matching 
(for symmetry constraint), the new element has much less degrees of freedom (dof) but a higher order of 
approximation property, compared to previous elements |2 H I30 [ l3T | . This is confirmed by numerical results. 
In 3D, the new element has only 15 stress plus 3 displacement dof on each element, much simpler than the 
first order element, with 54 plus 12 dof per element, of [23]. Notice that the element of [23] is previously 
the simplest rectangular element in 3D. 

The rest of the paper is organized as follows. The minimal element in 2D is introduced in Section 2. The 
well-posedness of the finite element problem, i.e. the discrete coerciveness and the discrete inf-sup condition, 
is proved in Section 3 for the pure displacement problem. The optimal order convergence is shown in Section 
4. The element is extended to any space-dimension in Section 5. In Section 6, the stability of the minimal 
element is shown for the pure traction problem. Numerical results in 2D and 3D, including that for a pure 
traction problem, are provided in Section [3 which show a superconvergence of the minimal elements herein. 

2. A MINIMAL ELEMENT IN 2D 

The 2D element is presented separately in this section for fixing the main idea while the whole family will 
be developed in Section 5. Also for simplicity we consider a pure displacement problem first. The analysis 
for other boundary value problems will be given in Section 6. 

Consider a pure displacement problem (and a pure traction problem in Section 6): 

(2.1a) div{A-\{u)) = f in n, 

(2.1b) M = on F = dn, 

The domain is assumed to be a rectangle (it is straightforward that results can be extended to domains 
which can be covered by rectangles) , which is subdivided by a family of rectangular grids Th (with grid size 
h). 

The set of all edges in Th is denoted by £h, which is divided into two sets, the set £h,H of horizontal edges 
and the set £h,v of vertical edges. Given any edge e £ £h^ one fixed unit normal vector n with components 
(ni, 712) is assigned. For each K € Th, define the affine invertible transformation 

Fk ■■ K ^ K, 

x\ fx 

y) \y 

with the center (sq.Kj Vo.k) of K, the horizontal length h^^K, and the vertical length hy^K, and the reference 
element i^ = [-1,1]^. 
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On each element K ^ Th, a, constant finite element space for the displacement is defined by 
(2.2) VK=Vo(K,R^) = {h) \vi,V2€Po{K)y, 

while the symmetric linear finite element space for the stress is defined by 

where subscript S indicates a symmetric matrix stress, and 

Pi.i{K) =span{l,a;}, 

■Pi{K) =span{l,a;, y}, 

A,2(^)=span{l,y}. 

The dimension of the space Vk is 2, and that of Sa' is 7. The nodal degrees of freedom for (wi,W2), ctu, 
and (722, are 

• the moment of degree on K for vi and V2 ; 

• the moments of degree on two vertical edges of K for an; 

• the moments of degree on two horizontal edges of K for (T22 ; 

The nodal degrees of freedom for ai2 will be studied as follows. Locally Vi{K) is the space of linear 
polynomials. Globally, let Wh be the Pi-nonconforming space on 7/i, which is first introduced in ^25^ as a 
nonconforming approximation space to H^{il.) on the quadrilateral mesh; see also [2D]. To be exact, W^ is 
the space of piecewise linear polynomials, which are continuous at all mid-edge points of triangulation Tn- 
Wh is the finite element space approximating function ai2. 

The global spaces T,h and Vh are defined by 

/CTli (712 \ ^ r2/ 



(2.4) Eh^{(J^( '' ''] € L'{n, §) I a\K e Ek for aU K e %, 

\(7i2 (722/ 

(7ii is continuous on all vertical interior edges, 
(722 is continuous on all horizontal interior edges, 
(712 is continuous at all mid-points of interior edges }, 

(2.5) Vh ^{v(E L^in^R^) I v\k e V{K) for all K e %}■ 

Since an is continuous on all vertical interior edges, the derivative dxan is well-defined in L-{Q). However, 
(712 is not continuous on 17 so that dxai2 and dyai2 are not in L-[Q). Therefore the discrete stress space 
Tih is a nonconforming approximation to iJ(div, f2, §). So the discrete divergence operator div/i is defined 
elementwise with respect to 7^, 

divft t\k = div(T|if ) Vt e S/^. 

The mixed variational form for (|2.1a|) is ()1.1|) . The mixed finite element approximation of Problem (|1.H) 
reads: Find [ah, Uh) G T,h x Vh such that 

f [Aah, t) + (div/iT, Uh) = Vt e S^, 

(2.6) < 

j^ (div,i ah,v) = (/, v) \/v eVh- 

It follows from the definition of E/^ that div^ Th are piecewise constant for any r^ G E/^, which leads to 

div/iSft C Vh- 
This, in turn, leads to a strong discrete divergence- free space: 

(2.7) Zh = {Th e ^h I (div;, Th,v)=Q yv(^Vh} 

= {t/i e S/i I div,i T/i = pointwise }. 
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For the analysis, define the following broken norm: 



(2.8) 



|div;,r||2)V2 VreE;,. 



iniff(divo = [\\T\i 

The rest of this section is devoted to an alternative definition to Wh, the space for (J12 in S^. The 
dimension of the space Vi {K) is three, less than the number of edges or vertexes of element K. The discrete 
shear stress cri2 is still defined by four vertex-value functionals, which are not linearly independent though. 
A constraint can be posed on those four functionals if one defines a functional set N on Vi{K), cf. ^51 
Lemma 2.1]. 

Here the idea from [20^ of a frame for Vi{K) will be used. To this end, define the frame for the space 
'Pi{k) = span{l, £, y) by 

1 — X — ii , 1 



1 



x-y 

4 

x + ij 



1 



-1.1 



X — y 

1 ' 

x + ij 



This frame is depicted in Figure [3l 



-1.1 



1 — X + y 



1-x-y 




1 + X + y 



1 + X — y 



Figure 3. Four nodal (frame/not basis) functions of Wh on K. 



An interpolation operator ni2, from H'^{Vl) (i.e., some continuous functions) to Wh is needed. The 
interpolation on K is defined as 

ni2(Ti2 = ai2{x^j^,y^j.)(j)-i-i + ai2{x^j^,y^j^)(j)i,^i 

+ (3-12(x3^^, 1/3^^)01,1 + Cri2{x^J., 2/4^^)0-1,1, 

where the four vertexes are numbered counterclock wise, 

i^2,k^y2.k) ^ (I'-l)' 
(^3,a:' %,/$:) = (!'!)' 
i^4,k^yA^k) = (-I'l)- 
In the same fashion, the interpolation II12 is defined on all K € Th by 

(2.9) Tli2cri2{x,y) = ai2{xi^K,yi,K)(t>-i-i{Fj^^{x,y)) 

+ (Ji2{x2,K,y2,K)(l)i-i{F^^{x,y)) 

+ cri2{x3^K,y3,K)(f>is{FK^{x, y)) 

+ ai2{x4^K,y4:.K)(f'-is{FK^ix,y)), 

where {x,y) £ K, and {xi^K,yi,K) are the four vertexes oi K. As (/)_i._i(0,— 1) = (/)i,_i(0, — 1) = 1/2, it 
follows that 

ni2cri2|e±(em) = - (cri2(ei) +0-12(62)) , 
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where e+ and e_ are two sides of an edge e € £h, Cm is the mid-point of e, and ei and ei are two endpoints 
of e. That is, ni2cri2 is continuous at all mid-points of edges. For a vertex in 7^, 

c,j = (iKjh), 0<t,j<N, N^ l/h, 

it may be shared by one, or two, or four elements K G Th- The combination of the frame functions at the 
vertex Qj- forms one global frame function 4>i,j. For example, at vertex co,i, as it is shared by two elements, 
Ki^i = [0, h] X [0, h] and Ki,2 = [0, h] x [h, 2h], 

r</.^i,i(f(x-|),f(y-|)) {x,y)€K,,,, 

V'oa - j 0-i,_i(f (x- |),f (y- f )) (a;,y) e Xi,2, 

lO elsewhere on ft. 

Note that ^/^ij- is not continuous at Cij. Thus, the finite element space for (T12 in ()2.4p is 

N 

(2.10) M^/, = {s e L2(f^) I s - 51 P^J^^.j}- 

3. Well-posedness of the discrete problem in 2D 

This section considers the well-posedness of the discrete problem (|2.6p . which needs the following two 
conditions. 

(1) K-ellipticity. There exists a constant C > 0, independent of the meshsize h such that 
(3.1) (AT,T)>C||r||2,(,i,,^) VrGZ,, 



where Zh is the divergence- free space defined in (|2.7p . 
(2) Discrete B-B condition. There exists a positive constant C > independent of the meshsize h, such 
that 

fon\ ■ f {diVhT,Vh) 

(3.2) ml sup T—j j— j— > C. 



Theorem 3.1. For the discrete problem i2.6]) . the K-ellipticity p.ip and the discrete B-B condition 
hold uniformly. Consequently, the discrete mixed problem (J2.6I) has a unique solution (cr/j, Uh) G S^ x V/j. 

Proof. It follows from (12. 7p that for all t e Z^, div/i r = 0. Thus |j div/iTJjo = and ||T||//(div^) = ||t||o. 
Since the operator A is symmetric and positive definite, the ii'-ellipticity of the bilinear form {At, t) follows. 
It remains to show the discrete B-B condition p.2p . Since the usual technique based on canonical 
interpolations operators for discrete stress spaces [3J[5] is inapplicable here, a constructive proof is adopted. 
For convenience, suppose that the domain fi is a unit square [0, 1]^ which is triangulated evenly into iV^ 
elements, {Kij}. For any v G Vh, it can be decomposed as a sum, 

N N 

(3-3) ^f>- = ^Yl ^'J "^^J (^' 2^) ' 

where (pij{x) is the characteristic function on the element Kij, and Vij = {Vi^ij, V2.ij) = {vh\Ki)- A discrete 
stress function r^ € S/j will be constructed with 

diYhTh = Vh and \\Th\\H{divh) ^ C'll^'/illo- 

The construction of Th is motivated by a simple proof of the inf-sup condition of the ID Raviart-Thomas 
element for the ID Poisson problem. The shear stress T12 can be taken zero, i.e., T12 = 0; the normal stress 
Til (resp. r22 ) of Th can be constructed so that it is independent of the second (resp. first) component of 
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Vh- In addition, th (resp. T22) can be a continuous piecewise linear function of the variable x (resp. y) and 
a piecewise constant function of y (resp. x). Therefore, they are of form 

■i-i 

(3.4) Tii{x,y) ^ h^ Vi,mj + Vi,ij{x - x^-i), 

771—1 

(3.5) T22{x,y) = h^V2,ik + V2,ij{y~yj-i), 

fe=i 

for Xi-i < X < Xi and yj-i <y<yj {{xi,yj) is the upper-right corner vertex of square Kij.) Thus, define 

By this construction, dxTn = (i^?i)i and dyT22 — {vh)2- This gives 

(3.6) div,,r/i = Vh- 
An elementary calculation gives 

N N 

\Vh\ 






N 

By the Schwarz inequality, 

kiillo = y2 I [ ^ y^ Vi^mj + Vi^ijix - Xi-i) I dxdy 



< 



j-^l-'^'J \ m=l 
JV „ / i-1 

J^l-^-li'J \ m=l 



i,mj) + iVi^tj){x - Xt-i) -idxdy. 



Further, since N = 1/h and /^., = h^, 

N / i \ N / N \ 

iiruiig <Y.[h"Y. (^1-^)' • ^/^^ < E ^' E (^1-^)' • ^''^^ 

i.j — 1\ 771—1 / j^l\ m— 1 / 

N 

= /^^E(^M.)'- 

A similar argument leads to 

llr22|l^</i'E(^2.,)'- 
The combination of the aforementioned two identities and two inequalities yields 

lk?j||ff(divO = \\Th\\l+ ||div,,T,J^ 

= ||rii||§ + ||r22||g + ||«„||^<2||t,„||2. 

Hence, for any Vh G Vh, the B-B condition (13. 2|) holds with C — 1/a/2: 

. (div/,T, w,0 . Il-y/illg 1 

ml sup -n— n n — tt- > mi — ;= ^— = —;=■ 

This completes the proof. 
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4. Error analysis in 2D 



The section is devoted to the error estimate stated in Theorem I4.3| which is based on the approximation 
error estimate of Theorem 14.11 and the consistency error estimate of Theorem 14.21 

In order to analyze the approximation error, for any r e _ff (div, Q, §) n H^{fl, §), define an interpolation 

fA -W TT /niiCTii ni2Cri2\ ^ 

(4.1) Uha^l e S;,, 

\lll2Cri2 li22Cr22/ 

where Tin and 1122 £ire standard, satisfying, respectively, 

(4.2) / IIiiCTiids = / (Tiids for any vertical edge e e £h, 

J e J e 

(4.3) / 1122 (722^5 = / o'22'^s for any horizontal edge e € £h- 

J e J e 

ni2 is the interpolation operator defined in (|2.9|) . from the space H^{n) to Wh- It is shown by Park and 
Sheen [35] that 

(4.4) \v-Ili2v\rn,K<Ch^-nvkK, m = 0,l, K ^ Th- 
Theorem 4.1. For any a € H^{fl,S), it holds that 

lk-n,,a|io<C/i|la|li, 



|divh(o- - n/^cr)||o < C/i||fT| 



2- 



Proof. By the scaling argument and the standard approximation theory, the following two estimates will be 
proved 

(4.5) \<Jii-nn<7ii\oM<Ch\an\i,K ^KeTh, 

(4.6) |_(aii-nnaii)|o,K<C/i|^|i,K "iK eTh- 

For any element K E Th, by (14.21) (i.e., the interpolation (14.21) is equivalent to a mid-point interpolation), 

llo-ii - Iliicriillo j^ = -^-^ / |(Tii - Iliianl^dxdy 
<Ch^\aii\lj^<Ch^\aii\lj,. 
This is (|4.5p . By the reference mapping, 



(4.7) 



Now 



d ' 

— (cTll -Iliicrii) 



^ -r- |TrT(o-ii -niio-ii)pdxdj/ 
hx Jk ox 

f d d - 



I -^tlii^iidxdy = / {(nii<Tii)(l,y) - (nii(Tii)(-l,2/)}d2/ 

= / {0-11(1,2/) -<7ii(-i,y)}dy 



-^7- dxdy, 
K ox 



SYMMETRIC MIXED FINITE ELEMENT 



This means -r— (HuO'ii) = P„ i>(—^ — ), where P^, ^ is the projection operator onto the constant space on 
ox ' ox ' 

element K. A substitution of it into (14.71) leads to 



d_ 
dx 



(fJii -Iliicrii) 



<C 



0,K 



d&i 



dx 






0,K 



<Cinf 

cGR 






By the Bramble-Hilbert Lemma, 



d 



— (cth -Iliicrii) 



<C 



0,K 



dan 



dx 



< Ch^ 



l.K 



0,K 



dan 



dx 



1,K 



This is gH). 

A similar argument yields 

(4.8) 
(4.9) 
Noting that the L^ norm on S is 



0-22 - n220-22||o,A' < C'/l|(T22 1 1,K VJST G Th, 



d 

T^(0'22 ^ 1122(722) 



< C/i 



O.K 



9(722 



5j/ 



VA' e %■ 



1,K 



Ikllo.K = Ikullo.K + 2||cti2|Io,a: + lk22|lo,A:, 
A combination of the estimates ()4.5p . (|4.6p . (|4.8p . ()4.9p and ()4.4p . completes the proof. I 

Theorem 4.2. Assume that {a,u) be the solution to the problem ()l.ip with u £ ifo(r2,M^) n iJ^(fi,M^). 

(A(7,Tft) + {<XwhTh,u) 



(4.10) 



sup 



< Ch\u\2. 



\\Th\\H{diY^) 

Proof. It follows from the first equation of (|l.ip that yl(7 = ^(Vu + Vu"^) for the exact solution u G 
i/Q(f7, M^). An elementwise integration by parts gives 

(e(?i),T/,) = -(divhT/i,ii) + V" / ThU-uds Vth E T,h, 

which implies 

(4.11) {Aa,Th) + {dWhTh,u) ^ y2 / Thu-uds. 

r^^^ JdK 



Ken 



-Til 

-T12 



ri2 
T22 



,e3,K 



msji 



e4,K *>mi,K X ""^^'^ 



— • — 



Til 
T12 



-T12 
-T22 



, ei.if 



, e2,_R- 



FiGURE 4. Th- n on the four edges of element AT, cf. (I4.12p . 



10 



JUN HU, HONGYING MAN, AND SHANGYOU ZHANG 



Let Tfi\K = { I J cf. Figure m Since rn is continuous in the x-direction and T22 is continuous 

.T21 T22, 



in the y-direction, there is a cancehation for these two components on the inter-element boundary. Since 

ueH^{n,R^)riH^{n, 



u,^A j, 



(4.12) 



y / T/i • nuds 

r^n- J OK 



K£Th 






)Ti2U2ds+{ - )Ti2Uids 

64, K Jei.K Je3,K 



For any v G H^{K), define the L^-projection operator Jg on an edge e by 

1 /■ 

JeV — Tl I vds. 

|e| Je 

Because T12 is continuous at the mid-point of all edges, it follows that, including boundary edges where 
Ui — 0, on the horizontal edges £h,H, 

^(/ -/ )Ti2Ulds = ^ j {Tl2\e+ ~ Ti2\e_)uids 

y^ j (Ti2\e+-Ti2\e_){ui- JeUi)ds. 



Ken 



eeEh,, 



After inserting a same constant JkTi2 — ^ r^ Ti2dxdy / \K\ into the two integrals on two horizontal edges of 
one element K, the sum can be rewritten as 



(4.13) 



^ ( / - / )Ti2Uids 

= ^ / Ti2(mi - Jei^Wl)ds- / Ti2{ui ~ Je^^Ui)ds 

= ^ / (ti2 - JkTi2){ui - Jei i^Ui)ds 



Ken- 



(ti2 - JkT12){ui - Je3^Ui)ds. 



There is some superconvergence property for the two terms in ()4.13p if they are considered together. In 
fact, on the reference element K, fi2(i, ±1) — fi2(0, 0) -I- xdxTi2{Q^ 0) ± 9t;fi2(0, 0), and Jj^ti2 — ti2(0, 0). 
The property of Je gives 



(t12 - JkTi2) (Wl - JeiUl){x, -1) - (ui - JesUi){x, 1) 



I f'd 



1 d 



1 



^"ic^y-o/ (""(^jl) -'"(i,-l))c^^ 






dx 
dx 
dx 



1 f'd 



X^T7TX2 



4y_i dx \-J_iJ_i dy 
AJ_^ dx lJ_j_J_j_Ji 



(-^"1(2:, 2/) - -^u{t,y)}dydt 
dy 



dx 



dxdy 



ui{s,y)dsdydt 



dx. 
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By the Schwarz inequality and (|4.13p . 
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^ ( / - / )Ti2Uids 



22 



Ken 



1 d 



2^ I pi^ 



Tl2X 



1 /■! r^ fp 



_i dx '-J-iJ-iJi dxdy 



:Ui{s,y)dsdydt 



dx 



\KeTh 



d_ 
dx 



Tl2 



Ch"- 



d_ 
dx 



Tl2 



n22 



o.kJ \Ken 

2\ 



52 



dxdy 



Ui 



0,K. 



52 



dxdy 



Ui 



< Ch^\n2\Uui\l. 
Here | • |i^;i is the elementwise seim-H^ norm. A similar argmxient bounds the other term in (I4.12p by 



^ ( / - / )Tl2U2ds 

KGTh ^2, ft: -^64, ft: 



<Ch^\Ti2\lM\u2\2- 



A combination of these two estimates with (14.111) implies 

\{Aa,Th) + {div hTh,u)\ < Ch^\u\2\Th\i.h- 
By the inverse inequality, 
(4.14) \{Aa,Th) + {diYhTh,u)\ < Ch\u\2\\Th\\o- 



Theorem 4.3. Let (a,u) Cz T, x V be the exact solution of problem p.ip and {Th,Uh) G S/j x V^ the finite 
element solution of (j2.6p . Then 

\\(J-(Jh\\o <C/l(||M||2 + ||cr||2), 

|ldiv,,(o- - o-,,)||o < Ch{\\u\\2 + Wah), 
\\u-Uh\\o<Ch{\\u\\2 + \\a\\2). 

Proof. Let 

Zf = {r e S,, I (div^r, v) = (/, v) ^v e V,,}. 

The finite element solution cr^ is in Zf. Thus, for any t € Zf, it holds Uh — t ^ Zh, i.e., 

div,i (cr/i -r) = 0. 

It follows from the ii:-ellipticity (cf. p.ip ) that, for all t E Zf, 

C\\(Jh -t\\1< {A{ah - T),(Th - r) 

= {A{a - t), ah~T) + {A{ah - a), ah ~ r) 

= {A{a -T),ah-T) - {Aa, (Th-t) - (div /^(cr/i ~T),Uh) 

= (A(cr - r), cr/i - r) - (Act, ct/^ - t) 

= (yl(CT - r), CT/i - r) - {Aa, ah - t) - {divhiah - t),u). 
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An application of the Schwarz inequality leads to 

\Wh -T\\H{d\Vh) = lk/l-T||o 

^ rn\ II (Act, (T/j -r) + (div /I (ct,, -r),-u) 

< <- |P - T||if (div;,) TTjT- -n 

== C'|lcr-r||ff(div„) + sup — — . 

TheSh L.|Fft||ff(divfc) 

By the triangle inequality, 

,.^K\ II II ^m- f II II , {Aa.Th) + {d\YhTh-u) 
(4-15) ||cr-fT,,||H(div;.) < C^i mf ||cr-r||H(div^)+ sup 1| 1| }. 

For a given Th € S;i, the discrete B-B condition p.2p ensures that the following problem has at least one 
solution 7;i g E;i, cf. [12], 

(4.16) (div,i7;i,w/i) = (div/i(cr - Th),Vh) V w/i G V^. 

It follows from the B-B condition (13.21) that 



1 (div,i7/i, Wft) 1 {diVhicr-Th, vh) 
7/1 H(divfc) < 7=; sup ^ — = — sup — 

< — ||div,,(cr-rft)||o. 
The identity (|4.16p asserts that •jh + Vi G Zf . The choice t = -fh + Th in (|4.15l) leads to 

\W -<^h\\H{divh) 

^nn\ II , II II , iAa,Th) + idivhTh,u) 

<C{\\a-Th\\H{divn) + \n\\H{divn)+ SUp . . } 

ThGSh Ir/i||ff(divh) 

^^M, II , {Aa,Th) + {divhrh,u) 
<C{\\a~Th\\H{divu)+ sup — }. 

TfcGSfe \\Th\\H(dWh) 

That is, 

iA T7\ II II ^ nf ■ r II II , {Aa,T) + {divhT,u) 

(4.17) \W-(Jh\\Hidiv,,) <C{ mf \\cr-Th\\H{dw,,)+ sup r— }. 

■^heSh ^gs,. \\T\\H{divh) 

The first term on the right-hand side of (|4.17p is the approximation error. The choice Th — Tlhcr with 
Theorem 14.11 gives its upper bound. The second term on the right-hand side of (J4.17I) is the usual consistency 
error for the nonconforming finite element method, which has already been bounded in Theorem 14.21 A 
combination of these two theorems implies 

\\(T-(Jh\\o <C/l(||M||2 + ||cr||2), 

\\divh{(r - CT/Ollo < Ch{\\u\\2 + hh)- 
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The rest of the proof is concerned with the estimation oi u — Uh- In view of the discrete B-B Condition 
p.2p . it holds, for any v €Vh^ 

C\\uh - v\\a 

(div,,T, Uh - v) (div,jT, Uh-u + u-v) 

- ^^P — \r~\\ = ^^P ^~n 

reSft Irllff(diVh) TGSh \\T\\H(diWh) 

[AwhT.Uh -u) 

< sup — — |-||u-i;||o 

TGEh \\T\\H{diVh) 

(divftT, Uh) + (Act, t) - (Act, t) - (div/jT, u) 

= sup T—r h \\U - v\\o 

reEh IrllH(diVft) 

(A(ct-ct/i),t) - (ACT,r) - (div,jT,u) 
= sup V\\u-v\\o 

rGEh |r||H(div,0 

(Act, r) + (div/iT, u) ^,,, ,, ,, ,, . 

< sup ^ '-^ ^^^ + C(||ct-ct,J1o + |1u-«|1o). 

By (|4.14p and the error estimation of ||ct — ct/,,||o, the triangle inequality plus v — Ph.u {Ph is the L^ projection 
into piecewise constant spaces) yield 

\\u - UhWo < \\u - Phu\\o + \\PhU - Uh\\o 

< Ch\u\2 + C(||ct -<Th\\o + \\u- PhuWo) 
<Ch{\\uh + \\ah). 

That completes the proof of this theorem. I 

5. The minimal element in any spatial dimension 

Assume the domain i7 is a unit hypercube [0, 1]" in the n-dimcnsional space, which is subdivided by a 
uniform rectangular grid of N"^ cubes: 

% := {^ii,i2,...,j,x = [(n ~ l)h,iih] X • •• [{in - l)h,i„h], 
l<iu---in<N; h=l/N}. 

The set of all {n — l)-dimensional face hyperplanes of the triangulation 7h that are perpendicular to the 
axis Xi is denoted by £n~i,i- That is 

£n-i,z = {[(«! - 1)/J, hh] X • • • X [{ii-i - l)h, ii^ih] x {i^h} 

X • • • X [{ir^ - l)h, inh], I < ii , ■ ' ' in < N , < ii < N} . 

The internal hyperplanes are denoted by 

£n-l,i{^) = £n~l,i H fl. 

The set of all (n — 2)-dimensional mid-surface hyperplanes (orthogonal to both Xi and Xj axes) are denoted 

by 

[{in-l)h,inh], l<ii,---^n<N, 0<H<N} 



£n~2,ij = {[{ii - l)h,iih] X ••• X {iih} x ••• x {{ij - -)h} x ••• x 



U {[(zi — l)h, iih] X • • • X {{ii )h} x • • • x {ijh} x • • • x 

[{in - l)h, inh], l<ii,---in<N, 0<1^<N}. 
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In addition, define £n-2,ij{K) := £n-2.ij n dK for any K ^Th- In 2D, these sets are 

1^1,1 = {all edges in Th perpendicular to xi }, 
iSi,2 = {all edges in Th perpendicular to X2 }, 
^0,12 = {all mid-points of edges in Th }■ 
In 3D, they are 

£2,i = {all squares in Th perpendicular to x^ }, 1 < J < 3, 
£i,ij = {all mid-square edges of squares in £2,1 and £2,j, 
parallel to x^ }, i y^ j y^ k £ {1, 2,3}. 
In n space-dimension, the symmetric tensor space is defined in (|1.2p . The discrete stress space is defined 

by 

(5.1) Tiilif G span{l,Xi},Tii is continuous on Ei e £n~i,i', 

TijIk & span{l,Xi,a;j}, nj is continuous on Eij e £n-2Aj{^) (■ 

Some comments are in order for this family of minimal finite element spaces. 

Remark 5.1. The normal stress tu is a constant on each {n— 1)- dimensional hyper-plane Ei e £n-i^i. In 
addition, for the case n = 1, Y^h is 

{tii G L (J7,R) t\k G span{l,x} is continuous at the nodes } CZ H (fi), 

the ID Raviart- Thomas space, which is the only conforming space in this family. 

Remark 5.2. The dimension of the space 

^h,ij := {nj e i^(fi,M) I Tij\K e span{l,Xi,a;j}, 

Tij is continuous on Eij € £n-2.ij{^)} 

is 

N^'-^iin + if - 1) = iV" + 2N"-^. 

see |25j for more details for 2D. 

Let us give the local basis for th and but a local frame (not basis) for r^ on an element K :~ Kij^^i2^,,,^i^ £ 
Th- Define, for (xi, . . . , Xn) G K, 



where 



V>l.^(xi,...,x„)=^W( w^^ / ^ )^ fc^o,l. 



^(o)(^)^i_i^ ^(i)(.i) = i±^, xe[-i,i]. 



2 

Define, for fc = 0, 1, 2, 3, for {xi, . . . ,x„) G K 



0!f^(.i,...,.„)^^W^ "''(^''^/'^^ "^ - ^'^ ' ^/'^^ 



where (cf. Figure[3]), for {x,y) G [—1,1 



V,- 

2 



h/2 ' h/2 
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Note that the above four functions are not linearly independent. In fact, 

0(0) _ 0(1) + 0(2) _ 0(3) ^ Q 

Then the finite element space can be alternatively defined by 

(5.2) S, = { (r,,)„,„ e L2(f7,M"X") | r,, = r,,; 

1 



fc=0 
3 



\K = Y.P^M%iK))cj^l^!Ki^U---,^n)}. 



fe=0 



p(fe) 



Here Tii{E^_^ ^{K)) are the values of t^ at the centers of the {n — l)-dimensional hyperplanes of K 
K ■ 

^^21 ,■■■ ,Zti ■ 



Ei'lUK) 



{ii-i - ^)h 
(ii ~ k)h 



fc = 0,l; 



Pij (-^1-2 (^)) € K are some parameters associated to the center-point of four (n— 2)-dimensional hyperplanes 
of K which are continuous on the four (two on the boundary) n-cubes sharing the point: 

fiH-k)h\ ({n-\)h\ ({ii~\)h\ ({H^\)h\ 



Ei%{K) 



i^^ - 0)h 



{H-l)h 






{i, - l)h 



V(^„-i)V \{ln~\)h) \{tn~^}hj \{t,,-l)hj 

As in 2D, the discrete displacement space is 

(5.3) Vfi ^ {v E _L"'(ri,M") I v\k is a constant vector }. 

In particular, the dof of the 3D mixed element is plotted in Figure [2j 

In the n-dimension, since div/i E/j C V/j, the K-ellipticity p.ip is proved exactly the same way as in 
2D. The explicit construction proof of the discrete B-B condition p.2p can be divided into n essentially 
1-dimensional construction proofs similar to that for the ID Raviart-Thomas element of the ID Poisson 
equation, see Section 3 for more details for 2D. For the consistency error in (|4.15p . the proof remains the 
same except there is a multiple summation instead of 2-index summation. All the analysis in 2D remains 
the same for n-D. 

6. The pure traction problem 

This section considers the pure traction problem, i.e., the stress space is subject to zero Neumann 
boundary condition while no boundary condition on the displacement. In practice, part of elasticity body 
should be located, i.e, the displacement has a Dirichlet boundary condition on some non-zero measure 
boundary. But the pure traction problem is the most difhcult one in mathematical analysis. A similar 
proof for Theorem 16.11 can prove it for partial displacement problems. For ease of presentation, details are 
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presented only for two dimensions. Note that the argument in any dimension is similar. The main idea is 
to use the macro-element technique where we construct a mass-preserving quasi-interpolation operator. 
Let RM be the rigid motion space in two dimensions, which reads 



RM := span ,„,,,.,,, 

1^ VOy \1) \—x 

Consider a pure traction problem: 

(6.1a) dYv{A-^e{u)) = / in 1] = (0, l)^ 

(6.1b) e{u) ■n = Q on T = dVi, 

(6.1c) (m,u) = VweRM. 

By the same discretization of uniform square grid Th with h ~ 1/N as in §2, the finite element equations 
(j2.6p remain the same except the spaces are changed with boundary and rigid-motion free conditions: 

(Act/i, r) -I- (div,i r, uh) = Vr e E/i,o, 

(div/i(Th,w) = (/, w) Vu e V/i,o, 

where 



(6.2) 



^h.o = {a=h^ ^1") e E;, I (7(me) • n - Vee(<f„nr)}, 



(6.3) 14,0 ^{v=r^\(^Vh\ {v,w) = Q Vwe RM}. 

Here me is the mid-point of an edge e, and S^ and Vh are defined in (12. 4p and ()2.5p . respectively. The 
earlier analysis remains the same except the discrete B-B condition p.2p as the stress space E/i,o is much 
smaller than before. 

Theorem 6.1. The following discrete B-B condition holds uniformly, 

mf sup - — — — > C. 

■"heVh.oa;.Gi;h,o Wh\\H(diY^)\\Vh\\o 

Proof. Let Vh = {{vh)i, {vh)2) = Y,i'"h)ijVi] as in ((3?3| . where [vhjij = {{vhji.ij, (i'/i)2,jj) is the constant 
value of Vh on square X^j. With the boundary condition on the stress, it is impossible to match {vh)i by 
dxTii alone as in (|3.4[) . That is, because the dof of {vh)i is n^ — 1.5 (due to a mixed constraint with the 
second component {vh)2), but the dof of {th} is only n{n — 1). This indicates that the help from dyTi2 is 
indispensable. But the traditional trick of interpolating smooth B-B stress function does not work here as 
(T12) does not have enough dof. In other words, the support of T21 is non-local, at least on four neighboring 
squares. Given Vh, a discrete B-B stress function will be constructed in two steps. First, a macro-element 
technique will produce a Uh globally so that Vh — div^ dh is rigid-motion free on each (2 x 2) macro-element 
K2i.2j,2h '■= [x2i,X2i+2] X [y2j,y2j+2]- In a sccond step, construct, macro-element by macro-element, a ah 
locally by internal dof only, so that divh o'h = Vh — div/j dh ■ 

To this end, define a local rigid-motion space on each macro-element K2i.2j.2h 

(6.4) R,j = spa.n{(t)l^^j,(j)'2.^j,(j)lij}, 

where 4>i ^ are defined in Figure[ni piecewise constant functions. Assume N is an even integer and decompose 
Vh into two parts, a local rigid-motion and a global rigid-motion-free part, 

(6.5) Vh = Vh + Vh, Vh ^ PL^(R,^)Vh ioT 0<i,j < N/2-1. 
Here the projection P£2(fl. .^w/i is defined as 

/ PLHRi^)Vh ■ if^^ijdxdy = / Vh ■ (f^^^jdxdy, to = 1, 2, 3. 
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iX2i,y2j) {X2i,y23) {X2i,y2j) 

Figure 5. Nodal values of three orthogonal basis functions, {0^ ^ -, m — 1, 2, 3}, of the 
rigid-motion space on macro-element K2i,2j,2h '■— [x2i,X2i+2] x [2/2^,2/2^+2], cf. (|6.4p . 

To construct (T;i, consider the pure traction PDE (16. lap with / — Vh with the solution u £ H^{Vt). Let 

Then 

(6.6) diver = Wh, ||cr||/f(div) < C'lI'D/illo. 
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Figure 6. Interpolation nodes for the Scott-Zhang C°-Qi Ih, ij^ , cf. (p?71) and /£, cf. (p^ . 

For the analysis, we need a mass-preserving quasi-interpolation operator. This will be achieved in four 
steps. First, let Ih be the boundary-condition preserving Scott-Zhang operator from [26;. which interpolates 
Hq functions to Cg(ri)-Qi(7/i) functions, shown in FigurelHl Then, we correct the mid-point values of edges 
of macro-elements to get a mass-preserving on each edge of each macro-element K2i,2j,2h- Let ttie be the 
mid-point of edge E of i^2i,2j,2h, which is also a vertex of Th, define the associated nodal basis function of 
the conforming bilinear element by 



^s("^b) = 1, Osiq) = for other vertexes q oiTh- 



Let 



ce ^ {v ~ Ihv)ds I / OEds. 
Define if : Hl{n) -> C^o{n)-Qi{Th) by 
(6.7) I^v = IhV + Y,CEeE- 

E 

Third, we correct the center value of I^v on each macro-element. Let ttLc be the center of macro-element 
K2i.2j.2h, which is also a vertex of 7^. Let the Qi nodal basis function 9ij for vertex nic be similarly defined 
as 6e- Define 



{v ~ l2jiv)ds + {v - l2JiV)ds) 

Ei,ij •'-E2,ij 



E2.ij 



,ds), 



18 JUN HU, HONGYING MAN, AND SHANGYOU ZHANG 

where Ei^ij — [0:21, 2^21+2] x {y2j+i} and E2AJ — {•'^21+1} x [y2j, 2/2^+2] are two intervals in the interior of 
K2i,2j,2h that take trie as their mid-points, cf. Figure [H Define /^ : Hq{Q) -^ Cg(rj)-(5i(7h) by 

(6.8) /;^z; = /fz; + ^c„%. 

FinaUy, define Hiz : i?o(^) ^ Wh,o := {w G Wh\ w{me) == Ve e f ^ n T} by 

(6.9) II12V := ni2/,> for any v € H^iQ). 

Since / Ili2l^vds — J I^vds for any e ^ £h , the definition of the interpolation operator ni2 leads to 



(6.10) 



/ IIi2vds = / vds and / (w — Ili2v)ds + {v — Ili2v)ds = 0, 

"E J E " E\^ij '' E2,ij 



for any E C dK2i^2j,2h and any macro-element i^2z,2j,2/i- In addition, 

(6.11) ||V,,ni2t^||o < C||Vt;||o. 
Then ct/j is defined as 

(6.12) CTii=niifTii, 

(6.13) 0-22 = n22Cr22, 

(6.14) 0-12 = ni2a'i2, 

where IIn, 1122 and ni2 are defined in (|4.2p . (14.31) and (|6.9p . respectively. 
We verify next, for ah defined in (|6.12p - (|6.14p . 



/ {divh ah - Vh) ■ (l)'in,ijdxdy = 0, to=1,2,3, 



for < ij < N/2. Note that div,, ct,, ^ Vh in general, though diver = Vh- From (li^ . (|i3|) . (|0|) and 
(|6.10p . and integrations by parts it follows 

X2i + 2 l'y2j + 2 

/ divh{a - ah) ■ (l)l^ijdydx 

X2i Jy2j 



Symmetrically, 



ry2j+2 

/ (I ~'nii)[aii{x2i+2,y) - (Jiiix2i,y)]dy 

■^y2j 

+ / (-/^ -ni2)[o-i2(a;, 2/2^+2) -0-12(2;, y2j)]rf2: = 0. 

J X2i 



a;2i+2 ry2j+2 

/ div;i (cr - cr,,) • (t>2^^jdy dx = 0. 

2=21 Jy2j 
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For the last preserved value, as diva = Vh pointwise, from (|4.2|) . (|4.3p . (|6.9p and ()6.10p . and integrations 
by parts it follows 

2^21 + 2 rV2j-i-2 

/ divh{cr - dh) ■ (j)l^ijdydx 

2/2 j + 2 

(/-nii)[crii(a;2i+2,y) ~ ox\{x2i,y)\dy 

y2j + l 

y2j+i 

(/ - nii)[(Tii(a;2i+2, y) - ox\{x2i,y)\ dy 
y2j 

X2i+1 

(I - Il22)[(^22{x,y2j+2) - 0-22(3;, y2j)] dx 

X2i 

a:2t+2 

(7-1122) [0-22 (a;, 2/2^+2) - (y22{x,y2j)]dx 

X2i+1 
X2i + 2 

{I -fli2)[ai2{x,y2j+2) + <7i2{x,y2j) - 2(Ji2{x,y2j+i)\dx 

X2i 
2/23 + 2 

(/- tii2)[2(Ti2{x2i+i,y) - (Ti2{x2i+2,y) - <J 12 {x2i , y)] dy = 0. 
y2j 

Thus 



Vh - divh CT/i 



K2 



J-R^j 0<i,j<N/2. 



We match next [vh — div/j ah] on each macro-element K2i^2j.2h by the divergence of internal 5 dof of discrete 

stress: 



'^ll,2i+l,2j+^^ '^ll,2i+l,2j+h '^12,2i+l,2j+^ ' ^22,2i+^,2j+l ^'^'^ ''': 



22,2i+f ,2j + l' 



where ^ii^2i+i,2j+^ denotes the value of cth at ((2i + l)/i, (2j + i)/i) and other notations are defined similarly. 
Note that the four mid-edge values of U12 are the same. Here on each macro-element, [vh — Arvh^h] is in 
the following space 



(6.15) 



M,j =span{0^,-, m = 1,2, 3, 4, 5} 



where 0^ j are defined in Figure [71 
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Figure 7. Nodal values of basis functions {0^^ ^ , 1 < m < 5} in M^j on macro-element 
X2i <x < a;2j+2, y2i < y < y2j+2, cf. (|6.15p . 

On each macro-element, define 5 stress functions to match the 5 basis functions of Af,-,- such that 



div;i (TmAj 
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Each such a function is denoted by a vector of its nodal values: 

' f''ll,2i+l,2j+l/2 

CTll,2i+l,2j"+3/2 
0'l2,2J+l/2,2j"+l/2 

^^■22,21+1/2,2^ + 1 
\ f''22,2i+3/2,2j + l / 



, ^m,ij 





(-l\ 




/o\ 




/-l\ 




/-l\ 




/o\ 









1 


1 


-1 


1 


-1 


1 





^=- 





7 





'2 


-1 


'2 


-1 


'2 

















-1 




-1 




-1 




\v 




\v 




\^) 




\-y 




lo/ 



for 1 < 771 < 5. A linear expansion 



Vh 



div/i dh 



Ko 



= E 



-'m,i]Vm.,ij 



defines ah{x,y) by 



<^h{x,y) 



5 

E 

7n— 1 



^ni^ij^ni^ij 7 v^^y) ^ -'^2i,2j?, 



2h- 



Thus 
(6.16) 



div,i ah ^Vh- divft ct^ and ||cr/i||o < C||t;/i||o. 



This stability is obtained by the standard scaling argument as all norms on 5-dimensional space Mij are 
equivalent. 

The final Uh for Vh is defined as 

<yh = CT/i + dh ■ 
As divft at = Vh, by ()6.1ip and (I6.16|) . the discrete B-B condition holds uniformly. I 

7. Numerical tests 

Two examples in 2D and one in 3D are presented to demonstrate the methods. These are pure displace- 
ment problem with a homogeneous boundary condition that u = on dQ. Assume the material is isotropic 
in the sense that 



(7.1) 



Aa = -— ( a- tT(a)S 

2fi\ 2/i + nA 



n = 2,3, 



1 0\ 
where ^ — I ^ 1 , and /i and A are the Lame constants such that < /ii < /i < /^2 and < A < cx). 



(7.2) 



(7.3) 



In 2D, let the exact solution on the unit square [0, 1]^ be 

4:x{l-x)y{l ~y) 
-4a;(l -x)y{l ~y)J ' 

and 

^e^~yx{l — x)y{l — y) 
sin(7ra;) sin(7ry) 

Notice that the second example is from [5T]. 
In 2D, the parameters in (|7.ip are chosen as 



A = 1 and 



1 

"=2- 



Then, the true stress function a and the load function / are defined by the equations in (jl.ll) , for the given 
solution u. 

In the computation, the level one grid is the given domain, a unit square or a unit cube. Each grid is 
refined into a half-size grid uniformly, to get a higher level grid, see the first column in Table [1] In Table [TJ 
the errors and the convergence order in various norms are listed for the true solution (j7.2p . Here and in rest 
tables in the section, Ih is the usual nodal interpolation operator. For example, IhUi{xi + h/2,yj + h/2) = 
Ui{Xt + h/2,yj + h/2), IhCTii{xi,yj + h/2) = aii{xi,yj + h/2), and Ihcri2 = ni2cri2, defined in (|2.9p . An 
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Table 1. The error and the order of convergence, for (|7.2I) . 
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Whu-UhWo 


/i" 


\\IhO' - <Jh\\o 


h" 


||div(4a-(7^)||o h" 


1 


0.05893 


0.0 


0.72887 


0.0 


1.41421356 0.0 


2 


0.02447 


1.3 


0.24585 


1.6 


0.35355339 2.0 


3 


0.00714 


1.8 


0.06587 


1.9 


0.08838835 2.0 


4 


0.00190 


1.9 


0.01708 


1.9 


0.02209709 2.0 


5 


0.00048 


2.0 


0.00440 


2.0 


0.00552427 2.0 


6 


0.00012 


2.0 


0.00113 


2.0 


0.00138106 2.0 


7 


0.00003 


2.0 


0.00029 


2.0 


0.00034526 2.0 



order 2 convergence is observed for both displacement and stress, see Table [TJ However, Theorem 14.31 onlv 
shows the first order convergence. Further studies on this superconvergence should be performed. 

The next example, (|7.3p . of Yi |31| is implemented for a comparison. The finite element errors and the 
order of convergence are listed in Table [2] An order 2 convergence is again observed. Notice that, see 
Figure [1] the minimal element of this paper has a much less dof than that of Yi, but has one order higher 
of convergence. 

Table 2. The error and the order of convergence, for (|7.3p . 





\\IhU-Uh\\o 


/i" 


\\IhO- - <Jh\\o 


h" 


||div(/^a-(7,,)||o h" 


1 


0.03619 


0.0 


3.08021 


0.0 


12.20143741 0.0 


2 


0.09843 


0.0 


0.54275 


2.5 


2.36338456 2.4 


3 


0.02594 


1.9 


0.15169 


1.8 


0.63139891 1.9 


4 


0.00664 


2.0 


0.03964 


1.9 


0.16050210 2.0 


5 


0.00167 


2.0 


0.01014 


2.0 


0.04029305 2.0 


6 


0.00042 


2.0 


0.00258 


2.0 


0.01008376 2.0 



As a third example, we compute a 3D solution for the following exact solution: 
(7.4) u 

on the unit cube [0, 1]^. This time, the parameters in (|7.1I) are taken as 



n6x{l ~ x)y{l - y)z{l ~ z)^ 

32x{l^x)y{l^y)z{l~z) 

v64x(l - x)y{l - y)z{l - z), 



A = 1, [I — — and n = 3. 
Again the order of convergence is still one higher than what is proved in this paper, see Table |3l 

Table 3. The error and convergence in 3D, for (|7.4p . 





\\IhU - u/i 


/i" 


//iCr — CT/i 


/i" 


||div(/,,a-(T^)||o /i" 


1 


0.16366 


0.0 


3.64496 


0.0 


8.94883415 0.0 


2 


0.07716 


1.1 


0.89446 


2.0 


1.73418255 2.4 


3 


0.02332 


1.7 


0.23153 


1.9 


0.42577123 2.0 


4 


0.00628 


1.9 


0.05946 


2.0 


0.10668050 2.0 


5 


0.00161 


2.0 


0.01518 


2.0 


0.02628774 2.0 
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As the last example, we compute the pure traction problem (|6.1ap with the exact solution 



(7.5) 



100x^(1 - xyy^{l - yy - - 



-1 



The matrix A is same as that in the first two examples. Our new finite element has no problem in solving 
the pure traction problems. The convergence results are listed in Tabled 

Table 4. The errors and the order of convergence for the pure traction problem (j7.5p . 





\\IhU-Uh\\o 


h" 


Who- - <Jh\\o 


h" 


|ldiv(/^a-a^)|!o h" 


2 


0.41470 


0.0 


1.19604 


0.0 


4.14320380 0.0 


3 


0.12546 


1.7 


0.26426 


2.2 


1.10584856 1.9 


4 


0.03273 


1.9 


0.06572 


2.0 


0.28799493 1.9 


5 


0.00827 


2.0 


0.01648 


2.0 


0.07297595 2.0 


6 


0.00207 


2.0 


0.00412 


2.0 


0.01830958 2.0 


7 


0.00052 


2.0 


0.00103 


2.0 


0.00458156 2.0 



[lo: 
[11 

[12 
[13 

[14 

[15 

[16 

[17] 

[18 

[19 

[2o; 
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